Image reconstruction device, image reconstruction method, image reconstruction program, and ct apparatus

ABSTRACT

A computerized tomography apparatus and program for obtaining a cross-sectional image corresponding to projections are provided in which, for a temporary cross-sectional image f(x, y) obtained in some manner, an evaluation function E is defined which includes differences between projections calculated from f(x, y) and measured projections, and f(x, y) is changed in a manner which substantially decreases E. The computerized tomography apparatus and program are characterized in which a back projection operation, which is required by conventional computerized tomography, is not essentially required. The computerized tomography apparatus and program are particularly effective in removal or reduction of metal artifacts, aliasing artifacts and the like.

TECHNICAL FIELD

The present invention relates to a technique of reconstructing across-sectional image of an object from the radiographic projections ofthe object.

BACKGROUND ART

Computerized tomography (CT) is a technique of reconstructing across-sectional image of an object from the radiographic projections ofthe object. FIG. 1 is a diagram schematically showing a typical X-ray CTapparatus.

In the typical CT apparatus, an X-ray source is moved around a targetobject, and irradiates X-rays to obtain the projections of the targetobject in many different directions. A cross-sectional image is obtainedby subjecting the projections thus obtained to a computationaloperation, so-called reconstruction. The Filtered Back-Projection method(FBP) is commonly used to reconstruct a cross-sectional image from theprojections. FBP is a kind of a transformation operation. In FBP, theprojections are subjected to a filtering essentially equivalent to thedifferential filtering, followed by “back projection,” in which eachprojection is projected back along the original projection direction,thereby a cross-sectional image is obtained. In this case, thedifferential filtering usually amplifies noise or errors, which can bethe source of artifacts (errors or false images which do not actuallyexist). Moreover, the back propagation operation spreads the artifactsthus created allover the cross-sectional image. Therefore, in CT, theartifacts often are not limited within a local portion around the sourceof the artifacts, and impairs the entire cross-sectional image,resulting in a fatal flaw.

Most artifacts are caused by the filtering operation and/or the backprojection operation involved in FBP. Therefore, if FBP is not used, across-sectional image substantially can be free from most of artifacts.As a method of calculating a cross-sectional image other than FBP, theAlgebraic Reconstruction Technique (ART) is historically important. ARTwas a major reconstruction method before FBP was proposed. In ART, theprocess of reconstruction is considered as a fitting problem where thecross-sectional image is a parameter and the projections are a targetdataset to be fit. The cross-sectional image is iteratively modified sothat projections (p) calculated from the cross-sectional image fitprojections (p₀) experimentally obtained. A feature of ART is that across-sectional image is asymptotically modified so that (p−p₀) becomeszero. ART usually requires a vast computation time in comparison to FBP.Therefore, ART is currently used only for particular applications (theanalysis of seismic waves, etc.). Although ART does not produce asextreme artifacts as FBP does, FBP often provides a more naturalcross-sectional image than ART does.

Besides the filtering operation and the back projection operation,artifacts may be caused by a lack or shortage of data in projections. Itis known that a lack or shortage of data often results in a fatalartifact especially in FBP. Other reconstruction techniques based onfitting, such as ART, are expected to be more robust against a lack orshortage of data than FBP. However, a lack of data is known to make CTan extremely “ill-posed problem,” under which it is essentiallydifficult to obtain reasonable solutions. One of the reasons why ARToften fails in fitting is that ART uses (p−p₀) for the target functionof fitting. So it is quite natural to consider the use of the (p−p₀)²instead of (p−p₀). In these cases, the least square method is one of themost popular way to minimize (p−p₀)². In the least square method, theinversion of a square matrix whose elements on one side is equal to thenumber of parameters is employed. Parameters in CT are values of pixelsin the cross-sectional image, and therefore, the number of theparameters is huge, If a cross-sectional image has 1000×1000 pixels, thenumber of the parameters becomes a million, and the number of elementsin the matrix is as huge as a trillion. Therefore, if the ordinary leastsquare method is used, the matrix is too huge to calculate. Instead ofthe ordinary least square method, the Simultaneous IterativeReconstruction Technique (SIRT) and the Iterative Least Square Technique(ILST) have been proposed. In these techniques, the calculation of across-sectional image is considered as a fitting problem as in ART. FBPis utilized as an inverse operation to calculate the cross-sectionalimage from projections partway through the calculation in both SIRT andWLST so as to circumvent the use of the ordinary least square method asdescribed above. Therefore, none of SIRT and ILST does not substantiallysolve the problems involved with the filtering operation and the backprojection operation as in FBP. This is probably why there have beenreports that SIRT and ILST just “reduce” artifacts.

Non-Patent Document 1: Yazdi M., Gingras L., Beaulieu L., An adaptiveapproach to metal artifact reduction in helical computed tomography forradiation therapy treatment planning: experimental and clinical studies,Int J Radiat Oncol Biol Phys, 62: 1224-1231, 2005.

DISCLOSURE OF THE INVENTION Problems to be Solved by the Invention

As can be seen from the discussion above, no use of FBP and robustnessagainst a lack of data are important so as to obtain a cross-sectionalimage free of artifacts. To achieve this, it is easily contemplated thata square error may be used as an evaluation function for fitting withoutusing FBP as in ART. Nevertheless, the artifact problem of CT has notbeen substantially solved, since this simple idea cannot bestraightforwardly realized. Firstly, the amount of calculation of across-sectional image by fitting is too huge to calculate quickly. It isessentially required to appropriately speed up the calculation.Secondly, the least square method, as a fitting algorithm, is weak (noteffective) against “ill-posed” problems. Therefore, any extension of theexisting algorithms (ILST, etc.) based on the least square method wouldnot be able to solve the artifact problem. Thirdly, other existingalgorithms as well as SIRT and ILST cannot completely eliminate thenecessity of FBP.

Solution to the Problems

A first novelty of the present invention is to employ SimulatedAnnealing (SA) as a fitting method. SA is a known technique, whichrequires a long time to perform fitting, and is inherently considerednot to be not suitable for CT. Despite this, by the sue of SA in CT, across-sectional image can be calculated by minimizing square errorswithout using FBP. SA is also stable even when fitting is performedunder ill-posed conditions, such as a lack of data or the like. Also inthis regard, the present invention is advantageous in terms of afundamental solution to artifacts. In view of the properties above, thepresent invention has an important novelty that a fundamental solutionto artifacts is obtained by applying SA to CT. When SA is simply appliedto CT, we need to iterate the calculation of projections from across-sectional image so many times. The amount of the calculation ofprojections from a cross-sectional image is substantially the same asthat of FBP. When SA is applied to CT, it takes several million times aslong a time as that of FBP. The calculation time can be of the order of“years” even if a state-of-the-art high-speed computer is used. In thepresent invention, this problem is solved by significantly reducing theamount of calculation by transformation of expressions. This solution isa technique required when SA is applied to CT, which is an importantfeature of the present invention.

A second novelty is the introduction of a smoothing term and an entropicterm, which actively destroy artifacts, into an evaluation function forfitting in addition to square errors. In conventional CT, only adifference between p and p₀ or their square errors is employed. In thiscase, there is always a possibility that a more satisfactory result(small errors) of fitting can be achieved while artifacts are left.Actually, this is often the case. Specifically, it is commonly expectedthat artifacts would be canceled with other artifacts. That is one ofthe reasons for the difficulties to eliminate artifacts. Even in thepresent invention, if only square errors are included in an evaluationfunction, artifacts were not completely eliminated. This fact indicatesthat another term is desired in addition to square errors so as toobtain a cross-sectional image free of artifacts. In the presentinvention, an entropic term and a smoothing term are introduced on thebasis of statistical mechanics. These terms mathematically represent anatural requirement that a cross-sectional image should be a “smooth”and “natural” gray image. The entropic term induces fitting so as todestroy artifacts and uniformize image quality of an entirecross-sectional image. The smoothing term suppresses the granularpattern of a cross-sectional image which is caused by the entropic term.By introducing the both terms, a natural cross-sectional image free ofartifacts can be obtained. Although the entropic term and the smoothingterm can each reduce artifacts singly, a combination thereof is found tobe more effective.

Note that the term “beam of radiation” as used herein refers toelectromagnetic waves, such as X-rays, visible light and radio waves,particle beams including electrons or charged particles, sound waves,which are vibration of a medium, and the like, in a broader sense thanthe general definition.

EFFECT OF THE INVENTION

A first effect of the present invention is to reduce artifacts which arecaused by a lack of data. Examples of a problem caused by a lack of dataincludes a case where an object to be observed has an opaque portion, acase where there is a limit on a projection angle, a case whereprojection angular intervals are not uniform, a case wherethree-dimensional CT (cone beam CT, helical scan CT, etc.) is performed,and the like.

In particular, it has been demonstrated that metal artifacts whichappear when an object to be observed has an opaque portion can bereduced. The term “metal artifact” means that when there is an opaqueportion (in many cases, a metal portion) with respect to X-rays in anobject to be observed, an entire cross-sectional image (not only theopaque portion) obtained by CT is destructively disturbed. Metalartifacts are caused by a discontinuous change in luminance of aprojection at an opaque portion, and a lack of information at the opaqueportion. When the differential filtering is applied to the opaqueportion in course of FBP, a discontinuous change in luminance takes anextraordinary value. Then, the extraordinary value is radially extendedto be a streak artifact via the back protection operation. Moreover, alack of information causes an unexpected contrast at the portions whichare not directly related to the opaque portion. Since the presentinvention does not use FBP and employ SA which is stable against thelack of information, it is easily understood that the present inventionis effective in removal of metal artifacts.

The present invention may be particularly useful for cone beam orhelical scan in terms of practical use. Both are calledthree-dimensional CT, and are currently rapidly becoming widespread.However, it is known that peculiar artifacts appear in three-dimensionalCT. The causes of the artifacts have been identified, but have not beensolved. In the case of cone beam, the cause of the artifacts is a lackof data. In the case of cone beam, conditions under which a completecross-sectional image can be obtained cannot be satisfied, so thatartifacts appear due to a lack of data. In the case of helical scan, thefundamental cause of artifacts is a back projection operation. Thegeometric anisotropy (helix) of the system of helical scan affects afiltering operation and a back projection operation, resulting inwindmill artifacts. Since the present invention does not require afiltering operation or a back projection operation, and is robustagainst a lack of data, the present invention can solve the problemswith three-dimensional CT.

Examples of a case where projection angles are not uniform includeanalysis of Earth's interior by CT using seismic waves, and analysis ofan atmospheric state by CT using radio waves from an artificialsatellite. These are known as typical cases where FBP cannot beutilized. It is expected that the present invention may be able toimprove the accuracy of analysis.

A second effect of the present invention is to increase a rate at whichprojections are taken and decrease the dose of X-rays. Since SA isstable against a lack or shortage of data, the present invention alsoinherits this feature. In the case of CT, the case where the number ofprojection angles is small is one of the cases of a shortage of data.Thus, the utilization of the present invention can reduce the number ofprojection angles as compared to conventional CT. The number ofprojection angles corresponds to the number of projections in which theobject is irradiated. So, the number of projections is proportional tothe imaging time and the dose. Therefore, a reduction in the number ofprojections leads to a decrease in the imaging time and the dose.

There is also a shortage of data when the image quality of projectionsis poor (low S/N ratio), for example. It has been demonstrated that thepresent invention is relatively stable even in such a case. If adecrease in the image quality of projections can be tolerated, this alsoleads to a reduction in the imaging time and the dose. Also, the presentinvention would contribute to an improvement in image quality in SPECTand PET in which the S/N ratio is extremely low.

A third effect of the present invention is that a luminance value of across-sectional image can be determined with high accuracy. The presentinvention provides a cross-sectional image which substantiallyfaithfully reproduces the measured projections. The accuracy ofreproduction is higher by about two orders of magnitude than that ofFBP. This is a benefit of the fitting algorithm which minimizes squareerrors. The higher accuracy of determination of a luminance valueguarantees the quantitativeness of the cross-sectional image and allowsa measurement of density using the luminance value. This feature can beused for an improvement in accuracy of measurement of bone density.Also, this would contribute to an improvement in accuracy of detectionof a pathological change (an organ containing a tumor, etc.).

A fourth effect of the present invention is that a cross-sectional imageobtained by the present invention has a higher contrast than that ofFBP. As described in the third effect, the present invention has thehigh accuracy of determination of a luminance value. As its subsidiaryeffect, the contrast of a cross-sectional image becomes higher. A highercontrast tends to lead to a higher apparent spatial resolution. As aresult, a cross-sectional image obtained by the present invention hashigher image quality than that of the conventional art. Notably, thiseffect allows the present invention to be useful for not only CT underspecial conditions or for special applications, but also ordinary CT.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic illustration of an X-ray CT apparatus.

FIG. 2 is a schematic diagram showing a relationship between across-sectional image f(x, y) and a projection p(r, θ).

FIG. 3 shows a typical example of p₀(r, θ), where the abscissa andordinate are respectively the projection angle and the channel portionof a detector.

FIG. 4 is a flowchart showing a basic procedure of the presentinvention.

FIG. 5( a) is a schematic diagram of a fitting process of steps (I) to(VI). FIG. 5( b) is a schematic diagram of a fitting process of steps(1) to (9).

FIG. 6 is the result of a technique described in Non-Patent Document 1.

FIG. 7 is a sinogram obtained from FIG. 6( a) by simulation (a whiteportion is an opaque region).

FIG. 8( a) is an enlarged diagram of FIG. 6( c). FIG. 8( b) is a diagramshowing a result of the present invention.

FIG. 9 is a schematic diagram showing an effect of each term in avirtual energy E.

FIG. 10 is a diagram showing artifacts and an artifact reducing effectof the present invention in the case with a limit on an angle.

BEST MODE FOR CARRYING OUT THE INVENTION First Embodiment

FIG. 1 shows a schematic diagram of an X-ray CT apparatus according toan embodiment of the present invention. The X-ray CT apparatus includesan imaging unit and a computer. The imaging unit includes a light sourcefor irradiating a target object (usually X-rays), and a detector forX-rays transmitted through the target object. The imaging unit obtainsprojections by emitting X-rays toward the target object in manydifferent directions. The computer includes a controller for the entireX-ray CT apparatus, and an image reconstruction processor for generatinga cross-sectional image of a region of interest of the target objectbased on X-ray projections obtained by the imaging unit. Note that theconfiguration of FIG. 1 is common to this and the following embodiments,but the image reconstruction processors of each embodiment performsdifferent process.

The image reconstruction processor of this embodiment employs SimulatedAnnealing (SA) as a fitting method for obtaining a cross-sectional imagefrom projections. Firstly, a framework of Simulated Annealing (SA) willbe described. SA is a fitting algorithm derived from a Monte Carlomethod, and is characterized in that the fitting is performed based onrandom numbers, and in that virtual energy and virtual temperature arehandled in the analogy to thermodynamics. SA is a known technique, whichis performed in steps (i) to (vi).

(i) A parameter is randomly chosen, and then, the parameter is changedbased on a random number (random number).

(ii) Evaluation is performed after the changing. A virtual energy E isconsidered as an evaluation function. In typical SA methods, E is takenas the sum of square errors. The change in E between before and afterthe changing is represented by ΔE (evaluation).

(iii) If the result of evaluation is improved (ΔE<0), the changing isaccepted (changing).

(iv) If the result of evaluation is worsened, the changing is acceptedin a probability of exp(−ΔE/T).

(v) The temperature T is decreased by a small amount.

(vi) The procedure above is repeated from (i).

In SA, since the changing is accepted according to Boltzmann statisticsas indicated in (iv) if the result of evaluation is worsened, thepossibility of escaping from a local minimum of ΔE is secured.Therefore, the method can reach a global solution without trapping in alocal minimum, and therefore, can stably function even under ill-posedconditions for fitting. Also, by gradually decreasing T, the methodgradually approaches a global solution (soft landing). A set of (i) to(v) is referred to as a single Monte Carlo step. In SA, the Monte Carlostep is infinitely repeated, the fitting process progresses. A timerequired for the calculation is obtained by a time required for oneMonte Carlo step x a required number of times of the Monte Carlo steps.The required number of times of the Monte Carlo steps is proportional tothe number of parameters (or the degree of freedom).

Next, this embodiment will be described according to the claims. If CTis considered as a fitting problem, a fitting parameter is across-sectional image (f(x, y)). Data to be fit is projections (p₀(r,θ)) where r indicates a channel position of a one-dimensional detectorused in the imaging unit, and θ indicates a projection angle. Adefinition of coordinates is shown in FIG. 3. p₀(r, θ) is a set of datawhich is the measuring projections while changing the angle θ, and canrepresent a two-dimensional image if it is plotted in a graph where thehorizontal axis represents r and the vertical axis represents θ. Suchdata is referred to as a sinogram. A typical sinogram is shown in FIG.3.

Roughly speaking, CT is a technique of obtaining a cross-sectional imagefrom a sinogram. In this embodiment, square errors between a temporarycross-sectional image f(x, y) and the measured projections p₀(r, θ) areused. In order to calculate the square errors, projections (p(r, θ)) arecalculated from f(x, y) by:

$\begin{matrix}{{p\left( {r,\theta} \right)} = {\sum\limits_{s}{f\left( {{{r\; \cos \; \theta} - {s\; \sin \; \theta}},{{r\; \sin \; \theta} + {s\; \cos \; \theta}}} \right)}}} & \left\lbrack {{Expression}\mspace{14mu} 101} \right\rbrack\end{matrix}$

Expression 101 represents the calculation of a projection of f(x, y)along a direction s by summing on s in FIG. 2. By using p(r, θ) thusobtained, the sum of square errors H can be calculated as follows:

$\begin{matrix}{H = {\sum\limits_{r,\theta}\left\{ {{p\left( {r,\theta} \right)} - {p_{0}\left( {r,\theta} \right)}} \right\}^{2}}} & \left\lbrack {{Expression}\mspace{14mu} 102} \right\rbrack\end{matrix}$

In typical fitting, Expression 102 is directly used as the virtualenergy E. This embodiment is characterized in that a smoothing term andan entropic term are introduced in addition to H. The virtual energy Eis defined as follows:

E=H−TS+cσ  [Expression 103]

where T represents a virtual temperature (temperature parameter), Srepresents an entropy, σ represents a standard deviation of pixelvalues, and c represents a coefficient which represents the strength ofthe smoothing term. In Expression 103, TS is an entropic term and cσ isa smoothing term. The definitions and calculation methods of S and Cwill be described below. In this embodiment, a cross-sectional image iscalculated in accordance with a procedure as shown in FIG. 4 using theseexpressions.

-   -   Step (a): An evaluation function (hereinafter referred to as an        “energy”) (E₀) is obtained with respect to a temporary        cross-sectional image f(x, y) which is prepared in some manner        (ST10 and ST20).    -   Step (b): (x₀, y₀) and Δμ are selected using random numbers, and        a portion of the cross-sectional image f(x, y) is modified        (ST30).    -   Step (c): An evaluation function E₁ is obtained with respect to        the modified cross-sectional image f(x, y) (ST40 and ST50).    -   Step (d): A differential (ΔE) between the energy (E₀) and the        energy (E₁) is obtained (ST60).    -   Step (e): It is determined whether or not the modification is        accepted, based on an acceptance function (typically Boltzmann        statistics) using the differential (ΔE) and the temperature        parameter (T) (ST70 to ST110).    -   Step (f): Control returns to step (a) (ST120).    -   Step (g): The value of the virtual temperature (T) is changed        every time the number of iterations of steps (a) to (f) reaches        a predetermined value (ST130).    -   Step (h): It is determined whether or not the result of        determination in step (e) satisfies predetermined stop        conditions. If the result of this determination is positive, the        process is ended. Here, it is assumed that the change is a        “success” if ST80 or ST100 is executed, and is a “failure” if        ST110 is executed. If the probability of success becomes lower        than 10% (this value is adjustable appropriately), the process        is ended. An estimated cross-sectional image at the end of the        process is considered as a cross-sectional image of the target        object, which may be displayed on a display of the computer or        may be recorded into a recording medium.

A series of operations from step (a) to step (h) corresponds to claim12. In this embodiment, these operations are performed by the imagereconstruction processor of FIG. 1.

These operations correspond to the basic steps (i) to (vi) of SA asfollows: step (b) corresponds to (i); steps (a), (c) and (d) correspondto (ii); step (e) includes (iii) and (iv); step (g) corresponds to (v);and step (f) corresponds to (vi). Therefore, in this embodiment, SA isfaithfully applied to CT.

Note that the determination of whether to end the process in step (h)may not be performed in the image reconstruction processor of FIG. 1.Estimated cross-sectional images may be successively displayed on adisplay of the computer, and a user who views the images may instructthe computer to end the process.

Second Embodiment

The virtual energy E indicated by Expression 103 is calculated bycalculating the sum of a series with respect to s, r and θ viaExpression 101 and Expression 102. A triple integral (the sum of aseries) is calculated, which takes a considerably long time. In otherwords, it takes a long time to calculate one Monte Carlo step. Inaddition, CT has a huge number of parameters. As a result, a totalcalculation time required for execution of SA is of the order of yearseven if a state-of-the-art computer is used.

Therefore, in this embodiment, instead of calculating E, only ΔE, whichis the difference when a change is made, is mainly calculated asfollows:

ΔE=ΔH+cΔσ−TΔS  [Expression 104]

Now, a change in a temporary cross-sectional image f(x, y) isrepresented as Δf(x, y). The change Δf(x, y) is a cross-sectional imagewhich has a value of Δμ only at a coordinate point (x₀, y₀) and zeroelsewhere. A projection Δp(r, θ) of Δf(x, y) can be calculated by thesame method as that of Expression 101. By using Δp(r, θ), ΔH can becalculated as follows:

$\begin{matrix}{{\Delta \; H} = {{\sum\limits_{r,\theta}\begin{Bmatrix}{{p\left( {r,\theta} \right)} +} \\{{\Delta \; p\left( {r,\theta} \right)} -} \\{p_{0}\left( {r,\theta} \right)}\end{Bmatrix}^{2}} - {\sum\limits_{r,\theta}\begin{Bmatrix}{{p\left( {r,\theta} \right)} -} \\{p_{0}\left( {r,\theta} \right)}\end{Bmatrix}^{2}}}} & \left\lbrack {{Expression}\mspace{14mu} 105} \right\rbrack\end{matrix}$

This expression is transformed as follows:

$\begin{matrix}{{\Delta \; H} = {\sum\limits_{r,\theta}\begin{Bmatrix}{{\Delta \; p\left( {r,\theta} \right)^{2}} + {2\Delta \; {p\left( {r,\theta} \right)}}} \\\left\lbrack {{p\left( {r,\theta} \right)} - {p_{0}\left( {r,\theta} \right)}} \right\rbrack\end{Bmatrix}}} & \left\lbrack {{Expression}\mspace{14mu} 106} \right\rbrack\end{matrix}$

Since Δf(x, y) has a value only at (x₀, y₀), Δp(r, θ) has a value of Δμonly at r(θ)=x₀ cos θ+y₀ sin θ and zero elsewhere. Therefore, theexpression within the braces { } in Expression 106 has values only atr(θ)=x₀ cos θ+y₀ sin θ. Therefore, the sum of a series does not need tobe calculated with respect to both r and θ, and Expression 106 can beexpressed as follows:

$\begin{matrix}{{\Delta \; H} = {\sum\limits_{\theta}\left\{ {{\Delta \; \mu^{2}} + {2\Delta \; {\mu \left\lbrack {{p\left( {{r(\theta)},\theta} \right)} - {p_{0}\left( {{r(\theta)},\theta} \right)}} \right\rbrack}}} \right\}}} & \left\lbrack {{Expression}\mspace{14mu} 107} \right\rbrack\end{matrix}$

Importantly, Expression 107 is the sum of a series only with respect toθ. As a result, the amount of calculation can be significantly reduced.Since p and p₀ are digital images, interpolation with respect to r(θ) isnecessary in advance to the summation in Expression 107. Therefore, theexpression within the braces { } in Expression 107 cannot be furtherexpanded. Despite this, if Expression 107 is expanded while admittingerrors, the following expression is obtained:

$\begin{matrix}{{\Delta \; H} = {{M\; \Delta \; \mu^{2}} + {2\; \Delta \; \mu {\sum\limits_{\theta}\left\{ {{p\left( {{r(\theta)},\theta} \right)} - {p_{0}\left( {{r(\theta)},\theta} \right)}} \right\}}}}} & \left\lbrack {{Expression}\mspace{14mu} 108} \right\rbrack\end{matrix}$

Note that M represents the number of projection angles. If Expression108 is used instead of Expression 107, the value of ΔH has an error ofabout 1%. However, Expression 108 can be calculated more quickly thanExpression 107, and therefore, is highly useful in the presentinvention.

Next, calculation of Δσ will be described. The standard deviation ofluminance values in an area around the coordinate point (x₀, y₀) isrepresented by σ. The area around (x₀, y₀) is assumed to be d×d pixelsaround (x₀, y₀) (in this example, d=5). The standard deviation of thesepixels can be obtained as follows:

σ=√{square root over (

f(x ₀ ,y ₀)

−

f(x ₀ ,y ₀)

)}{square root over (

f(x ₀ ,y ₀)

−

f(x ₀ ,y ₀)

)}  [Expression 109]

where:

$\begin{matrix}{{\langle{f\left( {x_{0},y_{0}} \right)}^{2}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{f\left( {{x_{0} + i},{y_{0} + j}} \right)}^{2}}}} & \left\lbrack {{Expression}\mspace{14mu} 110} \right\rbrack \\{{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{f\left( {{x_{0} + i},{y_{0} + j}} \right)}}}} & \left\lbrack {{Expression}\mspace{14mu} 111} \right\rbrack\end{matrix}$

By using Expressions 109 to 111, As can be calculated as follows:

$\begin{matrix}{{\Delta\sigma} = \sqrt{\begin{matrix}{{\langle{f\left( {x_{0},y_{0}} \right)^{2}}\rangle} - \frac{f_{i}^{2} - f_{j}^{2}}{d^{2}} -} \\\left\{ {{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} - \frac{f_{i} - f_{j}}{d^{2}}} \right\}^{2}\end{matrix}}} & \left\lbrack {{Expression}\mspace{14mu} 112} \right\rbrack\end{matrix}$

Note that f_(i) and f_(j) represent values of f(x₀, y₀) before and afterchanging.

The definition of the entropy S will be described before the followingcalculation of ΔS will be described. In general, images handled bycomputers are digital images. Not only the coordinates (x, y) of pixels,but also the values of the pixels have digital values. Therefore, it canbe assumed that a pixel is a quantum and that a pixel value is a quantumstate. In this sense, an image can be considered as a kind of quantumensemble. According to statistical mechanics, the entropy of such asystem is defined as follows:

$\begin{matrix}{S = {k\; \ln \frac{N!}{{N_{1}!}{N_{2}!}\mspace{14mu} \ldots \mspace{14mu} {N_{i}!}\mspace{14mu} \ldots \mspace{14mu} {N_{n}!}}}} & \left\lbrack {{Expression}\mspace{14mu} 113} \right\rbrack\end{matrix}$

In Expression 113, N represents the total number of pixels in an image,N_(i) represents the total number of pixels having a pixel value of i,which is a digital value, and k represents the Boltzmann constant intypical physics, but here, any value since the present invention is notdirectly related to physics. In the present invention, S is also assumedto be defined in the area of d×d pixels around (x₀, y₀) as is similar toσ.

A differential ΔS of S defined in Expression 113 is considered. It isassumed that by changing, a pixel value (digital value) is changed fromi to j. In this case, ΔS can be written as follows:

$\begin{matrix}{{\Delta \; S} = {{k\; \ln \frac{N!}{{N_{1}!}{N_{2}!}\mspace{14mu} \ldots \mspace{14mu} {\left( {N_{i} - 1} \right)!}\mspace{14mu} \ldots \mspace{14mu} {\left( {N_{j} + 1} \right)!}\mspace{14mu} \ldots \mspace{14mu} {N_{n}!}}} - {k\; \ln \frac{N!}{{N_{1}!}{N_{2}!}\mspace{14mu} \ldots \mspace{14mu} {N_{i}!}\mspace{14mu} \ldots \mspace{14mu} {N_{j}!}\mspace{14mu} \ldots \mspace{14mu} {N_{n}!}}}}} & \left\lbrack {{Expression}\mspace{14mu} 114} \right\rbrack\end{matrix}$

Thus, by changing, N_(i) is decreased by one while N_(j) is increased byone. If Expression 114 is expanded, a considerably simple expression isobtained as follows:

ΔS=k ln N _(i) −k ln(N _(j)+1)  [Expression 115]

In summary, the procedure of this embodiment is as follows.

(I) Using random numbers, (x₀, y₀) and Δμ are selected.

(II) Using Expressions 104, 108 (or 107), 112 and 115, ΔE is calculated.

(III) If the result of evaluation is improved (ΔE<0), Δμ is added tof(x₀, y₀).

(IV) If ΔE>0, Δμ is also added to f(x₀, y₀) with a probability ofexp(−ΔE/T).

(V) The temperature T is decreased by a small amount.

(VI) The procedure above is repeated from (I).

(VII) If the predetermined stop conditions are satisfied, the process isended. For example, the stop condition may be that the probability ofsuccess becomes lower than a predetermined value (e.g., 10%, this valueis adjustable appropriately) where “success” means the case when Δμ isadded in (III) and (IV). An estimated cross-sectional image at the endof the process is considered as a cross-sectional image of the targetobject, which may be displayed on a display of the computer or may berecorded into a recording medium.

It is seen that steps (I) to (VI) faithfully correspond to basic steps(i) to (vi) of SA.

Note that the processes in steps (I) to (VII) are performed in the imagereconstruction processor of FIG. 1.

Note that the end of the process in step (VII) may not be performed inthe image reconstruction processor of FIG. 1. Estimated cross-sectionalimages may be successively displayed on a display of the computer, and auser who views the images may stop the process.

Third Embodiment

Next, an important variation of the aforementioned method will bedescribed. In the aforementioned method, Expression 107 or 108 is usedto obtain the sum of a series with respect to θ. Therefore, each MonteCarlo step includes M times of summation. On the other hand, thisembodiment is provided to further reduce the amount of calculation.

Firstly, a back projection g(x, y) of p(r, θ) will be considered.

$\begin{matrix}{{g\left( {x,y} \right)} = {\sum\limits_{\theta}{p\left( {{{x\; \cos \; \theta} + {y\; \sin \; \theta}},\theta} \right)}}} & \left\lbrack {{Expression}\mspace{14mu} 116} \right\rbrack\end{matrix}$

In Expression 116, since a filtering operation is not performed, g(x, y)becomes a similar image to the cross-sectional image f(x, y), butsignificantly blurred. If g(x, y) is used, Expression 108 is written asfollows:

ΔH=MΔμ ²+2Δμ{g(x ₀ ,y ₀)−g ₀(x ₀ ,y ₀)}  [Expression 117]

Note that g₀(x, y) is a back projection of p₀(r, θ) and is calculated inthe manner similar to that of Expression 116. Expression 117 is superiorto Expression 108 in terms of the absence of summation. Since g₀(x, y)does not change at all during the calculation process, g₀(x, y) can becalculated in advance. On the other hand, g(x, y) changes a little forevery time a point of f(x, y) is changed, and therefore, to be exact,needs to be recalculated every time when a set of steps (I) to (IV) isexecuted. However, if it is assumed that the change in g(x, y) due tothe change in f(x, y) is small, another set of steps is applicable,which is of this embodiment. The steps of this embodiment will behereinafter described. Note that steps (1) to (9) below are performed inthe image reconstruction processor of FIG. 1.

(1) g₀(x, y) is obtained from p₀(r, θ).

(2) p(r, θ) is calculated from f(x, y), and then g(x, y) is obtainedfrom p(r, θ).

(3) An image Δu(x, y) whose pixel value is a change value of f(x, y) isgenerated using a random number.

(4) Expression 117 is applied to each pixel value to calculate an imageΔH(x, y).

(5) Similarly, Expressions 112 and 115 are applied to each pixel valueto calculate images Δσ(x, y) and ΔS(x, y).

(6) ΔE(x, y) is calculated based on Expression 104.

(7) Δμ(x, y) is set to 0 for a coordinate point (x, y) having a positiveΔE.

(8) Δu(x, y) is added to f(x, y).

(9) T is multiplied by α (α<1), and the procedure is repeated from (2).

(10) If the predetermined stop conditions are satisfied, the process isended. For example, the stop condition may be that the probability ofsuccess becomes lower than a predetermined value (e.g., 10%, this valueis adjustable appropriately) where “failure” means the case when Δμ(x,y) is set to 0. An estimated cross-sectional image at the end of theprocess is considered as a cross-sectional image of the target object,which may be displayed on a display of the computer or may be recordedinto a recording medium.

Note that the determination of whether to end the process in step (10)may not be performed in the image reconstruction processor of FIG. 1.Estimated cross-sectional images may be successively displayed on adisplay of the computer, and a user who views the images may instructthe computer to end the process.

Schematic concepts of steps (I) to (VII) and steps (1) to (10) areschematically shown in FIG. 5. In FIG. 5, an axis represents a value ofa parameter (in this case, a value of a pixel), and the virtual energy Eis represented by a gray level (a higher gray level indicates a smallerE). For the sake of convenience, the number of parameters (the number ofpixels) is assumed to be two. In the case of steps (I) to (VII) shown inFIG. 5( a), a minimum value of E is eventually reached via a zigzag pathwhile sometimes directs a wrong direction. On the other hand, in thecase of steps (1) to (10) of FIG. 5( b), a minimum value of E issearched for while the gradient of E is taken into consideration. Themethod of FIG. 5( b) would be intuitively recognized as being moreefficient. This embodiment is advantageous in terms of calculation speedbecause of a reduction in the amount of calculation by Expression 117and an improvement in search efficiency shown in FIG. 5( b). Note thatsteps (1) to (10) do not faithfully correspond to the basic algorithm(i) to (vi) of SA, and are a natural expansion of steps (I) to (VII).

Other Embodiments

An image reconstruction processor for executing the process described ineach embodiment above can be implemented using a program for causing acomputer to execute these processes, a computer in which the program isinstalled, a specialized LSI for executing these processes, or the like.

EXAMPLES

As an example, the effect of removing a metal artifact using simulationwill be described. Firstly, for comparison, results from the techniqueof Non-Patent Document 1 are shown in FIG. 6 (FIG. 5 in Non-PatentDocument 1).

In FIG. 6, results (a), (d) and (g) in the left column are originalphantom images without a metal artifact. Results (b), (e) and (h) in themiddle column are images reconstructed by a typical FBP method, whereinvisible regions are set at predetermined positions in the originalphantom images. Results (c), (f) and (i) in the right column are imagesreconstructed by the method of Non-Patent Document 1, where opaqueregions are set as in results (b), (e) and (h). The results (c), (f) and(i) are one of the best achievements in removal of metal artifacts bythe conventional arts.

The image of FIG. 6( a) in Non-Patent Document 1 was used, and theopaque regions were set at the same positions as those in the document.Then, projections were calculated by simulation. The result is shown inFIG. 7. FIG. 7 corresponds to p₀(r, θ) of this example. An imagereconstructed from FIG. 7 by the present invention and an enlargeddiagram of FIG. 6( c) for comparison are shown in FIG. 8. The effect ofthe present invention can be clearly seen from FIG. 8. In FIG. 8( b) ofthis embodiment, it is difficult to find even a feature of metalartifacts.

Note that, in FIG. 8( b) of this embodiment, edge portions are slightlyblurred. This is because a “factor which smoothes a cross-sectionalimage” (smoothing term cσ) is strong to some extent. It is known thatthe factor needs to be set to be strong to some extent so as to reducemetal artifacts. In any case, according to this embodiment, metalartifacts can be considerably removed, although there is still room foran improvement in the balance between each factor when the virtualenergy E is calculated.

Note that, for reference, a result obtained by executing the algorithmof this embodiment without setting the smoothing term cσ or the entropicterm TS is shown in FIG. 9. FIG. 9( a) shows a result obtained in theabsence of the smoothing term cσ (in the presence of only the entropicterm TS), FIG. 9( b) shows a result obtained in the absence of theentropic term TS (in the presence of only the smoothing term cσ), andFIG. 9( c) shows a result obtained in the absence of both the entropicterm TS and the smoothing term cσ. When FIG. 9 is compared to FIG. 8(b), it is found that, in this embodiment, metal artifacts can be removedto a greater extent when both smoothing term cσ and the entropic term TSare introduced to the energy E than when only either of them is set andwhen none of them is set. When only the smoothing term is set, streaksextending radially from an opaque region are left as artifacts. On theother hand, when only the entropic term is set, a cross-sectional imagebecomes granular and has a low S/N. For final image quality, the ratioof the coefficient c of the smoothing term and the virtual temperature Tduring a slow cooling process seems to be important.

Next, an example with the limitation on a rotational angle is shown inFIG. 10. FIG. 10( a) shows a result from reconstruction only by FBP(typical CT). FIG. 10( b) shows a result from application of ILST (animage reconstruction method based on the least square method). FIG. 10(c) shows a result from reconstruction by this embodiment. When theseresults are compared, it is found that this embodiment is effectiveagainst artifacts which appear in the case of the limitation. Artifactsdue to the angle limitation are characterized in that a circular regionis altered to an almond-shaped region having two opposite sharp ends andluminance is reversed across round side edges connecting the sharp endsof the almond-shaped region. In this embodiment, both of thecharacteristics are reduced.

INDUSTRIAL APPLICABILITY

The present invention is significantly effective against metalartifacts, and therefore, is particularly useful in a field where metalartifacts are serious, such as CT for teeth, CT for an object includinga metal implant, or the like.

Also, the present invention is generally effective in reconstructionfrom a set of projections having a lack of information. For example,there is a significant lack of information when there is a limit on aprojection angle. The projection angle limitation causes a problem withthree-dimensional electron microscopy, CT mammography, translation CT(Japanese Unexamined Patent Application Publication No. 2006-71472) orthe like. The lack-of-information problem also occurs inthree-dimensional CT, such as cone beam CT, helical scan CT or the like.The present invention is also effective in removal of artifactsappearing in three-dimensional CT.

Moreover, the present invention is applicable to a system in which theamount of information is considerably small. For example, the presentinvention is useful for fluorescent X-ray CT, seismic CT for imagingEarth's interior, and the like.

The present invention also has a subsidiary benefit that a luminancevalue (corresponding to an absorption coefficient in the case of X-ray)in a cross-sectional image can be determined with higher accuracy thanthat of the conventional art, for example. This effect can be applied soas to improve the accuracy of measurement of bone density or the like.

The present invention can be used to obtain a reconstructed image havinga higher contrast than that of the conventional art. Therefore, thepresent invention is also highly useful for typical X-ray CT in whichartifacts or the like do not cause a problem. Also, since the presentinvention is stable even when there is a shortage of data, the presentinvention is effective in a reduction in time required for measurementof projections, and therefore, a reduction in X-ray dose. According tothese features, the present invention has the potential to replace allexisting X-ray CT techniques.

1-14. (canceled)
 15. An image reconstructing device for obtaining across-sectional image of an object from projections (hereinafterreferred to as “radiographic projections”) obtained by irradiating theobject with a beam of radiation, comprising: means (a) for obtaining anevaluation function (hereinafter referred to as an “energy”) (E₀)including differences between projections calculated from a currentestimated cross-sectional image of the object and the radiographicprojections; means (b) for modifying a portion of the current estimatedcross-sectional image; means (c) for obtaining an energy (E₁) includingdifferences between projections calculated from the modified estimatedcross-sectional image and the radiographic projections; means (d) forobtaining a differential (ΔE) between the energy (E₀) and the energy(E₁); means (e) for determining whether or not the modification is to beaccepted, based on an acceptance function using the differential (ΔE)and a temperature parameter (T) for controlling an acceptanceprobability, and reflecting a result of the determination on the currentestimated cross-sectional image; and means (f) for changing a value ofthe temperature parameter (T) every time the number of iterations of aseries of processes of the means (a) to (e) reaches a predeterminedvalue.
 16. The image reconstructing device of claim 15, wherein themeans (a) calculates an energy (E₀) including a sum of differencesbetween projections calculated from the current estimatedcross-sectional image and the radiographic projections, and a standarddeviation of a local region of the current estimated cross-sectionalimage, and the means (c) calculates an energy (E₁) including a sum ofdifferences between projections calculated from the modified estimatedcross-sectional image and the radiographic projections, and a standarddeviation of a local region of the modified estimated cross-sectionalimage.
 17. The image reconstructing device of claim 15, wherein themeans (a) calculates an energy (E₀) including a sum of differencesbetween projections calculated from the current estimatedcross-sectional image and the radiographic projections, and an entropyof a local region of the current estimated cross-sectional image, andthe means (c) calculates an energy (E₁) including a sum of differencesbetween projections calculated from the modified estimatedcross-sectional image and the radiographic projections, and an entropyof a local region of the modified estimated cross-sectional image. 18.The image reconstructing device of claim 15, wherein the means (a)calculates an energy (E₀) including a sum of differences betweenprojections calculated from the current estimated cross-sectional imageand the radiographic projections, a standard deviation of a local regionof the current estimated cross-sectional image, and an entropy of thelocal region of the current estimated cross-sectional image, and themeans (c) calculates an energy (E₁) including a sum of differencesbetween projections calculated from the modified estimatedcross-sectional image and the radiographic projections, a standarddeviation of a local region of the modified estimated cross-sectionalimage, and an entropy of the local region of the modified estimatedcross-sectional image.
 19. The image reconstructing device of claim 15,comprising: instead of the means (a), (c) and (d), means (h) forcalculating ΔH using [Expression 1], and obtaining ΔE (ΔE=ΔH+ . . . )including the calculated ΔH as a component, $\begin{matrix}{{\Delta \; H} = {\sum\limits_{\theta}\left\{ {{\Delta \; \mu^{2}} + {2\; \Delta \; {\mu \left\lbrack {{p\left( {{r(\theta)},\theta} \right)} - {p_{0}\left( {{r(\theta)},\theta} \right)}} \right\rbrack}}} \right\}}} & \left\lbrack {{Expression}\mspace{14mu} 1} \right\rbrack\end{matrix}$ where, when the current estimated cross-sectional image ofthe object is represented by f(x, y) and the portion modified by themeans (b) is represented by Δf(x, y), Δf(x, y) is a cross-sectionalimage having a value of Δμ only at a coordinate point (x₀, y₀) and zeroelsewhere, and p(r, θ) represents a projection calculated from thecurrent estimated cross-sectional image of the object, p₀(r, θ)represents a radiographic projection of the object, r represents achannel position of a one-dimensional detector taking the projection, θrepresents a projection angle, and r(θ)=x₀ cos θ+y₀ sin θ.
 20. The imagereconstructing device of claim 15, comprising: instead of the means (a),(c) and (d), means (h) for calculating ΔH using [Expression 1], andobtaining ΔE (ΔE=ΔH+ . . . ) including the calculated ΔH as a component,$\begin{matrix}{{\Delta \; H} = {{M\; \Delta \; \mu^{2}} + {2\; \Delta \; \mu {\sum\limits_{\theta}\left\{ {{p\left( {{r(\theta)},\theta} \right)} - {p_{0}\left( {{r(\theta)},\theta} \right)}} \right\}}}}} & \left\lbrack {{Expression}\mspace{14mu} 2} \right\rbrack\end{matrix}$ where, when the current estimated cross-sectional image ofthe object is represented by f(x, y) and the portion modified by themeans (b) is represented by Δf(x, y), Δf(x, y) is a cross-sectionalimage having a value of Δμ only at a coordinate point (x₀, y₀) and zeroelsewhere, and p(r, θ) represents a projection calculated from thecurrent estimated cross-sectional image of the object, p₀(r, θ)represents a radiographic projection of the object, r represents achannel position of a one-dimensional detector taking the projection, θrepresents a projection angle, r(θ)=x₀ cos θ+y₀ sin θ, and M representsthe number of projection angles.
 21. The image reconstructing device ofclaim 19, wherein the means (h) calculates Δσ using [Expression 3], andobtains ΔE (ΔE=ΔH+cΔσ+ . . . ) including as a component a sum of aproduct (cΔσ) of the calculated Δσ and a coefficient c, and the ΔH,$\begin{matrix}{{\Delta \; \sigma} = \sqrt{\begin{matrix}{{\langle{f\left( {x_{0},y_{0}} \right)^{2}}\rangle} - \frac{f_{i}^{2} - f_{j}^{2}}{d^{2}} -} \\\left\{ {{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} - \frac{f_{i} - f_{j}}{d^{2}}} \right\}^{2}\end{matrix}}} & \left\lbrack {{Expression}\mspace{14mu} 3} \right\rbrack\end{matrix}$ where σ represents a standard deviation of luminancevalues of d×d pixels around the coordinate point (x₀, y₀) and iscalculated by [Expression 4], and f_(i) and f_(j) represent values off(x₀, y₀) before and after the modification by the means (b),$\begin{matrix}{{\sigma = \sqrt{{\langle{f\left( {x_{0},y_{0}} \right)}^{2}\rangle} - {\langle{f\left( {x_{0},y_{0}} \right)}\rangle}^{2}}}{where}} & \left\lbrack {{Expression}\mspace{14mu} 4} \right\rbrack \\{{\langle{f\left( {x_{0},y_{0}} \right)}^{2}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{f\left( {{x_{0} + i},{y_{0} + j}} \right)}^{2}}}} & \left\lbrack {{Expression}\mspace{14mu} 5} \right\rbrack \\{{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{{f\left( {{x_{0} + i},{y_{0} + j}} \right)}.}}}} & \left\lbrack {{Expression}\mspace{14mu} 6} \right\rbrack\end{matrix}$
 22. The image reconstructing device of claim 20, whereinthe means (h) calculates Δσ using [Expression 3], and obtains ΔE(ΔE=ΔH+cΔσ+ . . . ) including as a component a sum of a product (cΔσ) ofthe calculated Δσ and a coefficient c, and the ΔH, $\begin{matrix}{{\Delta \; \sigma} = \sqrt{\begin{matrix}{{\langle{f\left( {x_{0},y_{0}} \right)^{2}}\rangle} - \frac{f_{i}^{2} - f_{j}^{2}}{d^{2}} -} \\\left\{ {{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} - \frac{f_{i} - f_{j}}{d^{2}}} \right\}^{2}\end{matrix}}} & \left\lbrack {{Expression}\mspace{14mu} 3} \right\rbrack\end{matrix}$ where σ represents a standard deviation of luminancevalues of d×d pixels around the coordinate point (x₀, y₀) and iscalculated by [Expression 4], and f_(i) and f_(j) represent values off(x₀, y₀) before and after the modification by the means (b),$\begin{matrix}{{\sigma = \sqrt{{\langle{f\left( {x_{0},y_{0}} \right)}^{2}\rangle} - {\langle{f\left( {x_{0},y_{0}} \right)}\rangle}^{2}}}{where}} & \left\lbrack {{Expression}\mspace{14mu} 4} \right\rbrack \\{{\langle{f\left( {x_{0},y_{0}} \right)}^{2}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{f\left( {{x_{0} + i},{y_{0} + j}} \right)}^{2}}}} & \left\lbrack {{Expression}\mspace{14mu} 5} \right\rbrack \\{{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{{f\left( {{x_{0} + i},{y_{0} + j}} \right)}.}}}} & \left\lbrack {{Expression}\mspace{14mu} 6} \right\rbrack\end{matrix}$
 23. The image reconstructing device of claim 19, whereinthe means (h) calculates ΔS using [Expression 7], and obtains ΔE(ΔE=ΔH−TΔS+ . . . ) including as a component a sum of a product (−TΔS)of the calculated ΔS and the temperature parameter (T), and the ΔH,ΔS=k ln N _(i) −k ln(N _(j)+1)  [Expression 7] where S represents anentropy of a local region image of d×d pixels around the coordinatepoint (x₀, y₀) and is calculated by [Expression 8], $\begin{matrix}{S = {k\; \ln \frac{N!}{{N_{1}!}{N_{2}!}\mspace{14mu} \ldots \mspace{14mu} {N_{i}!}\mspace{20mu} \ldots \mspace{14mu} {N_{n}!}}}} & \left\lbrack {{Expression}\mspace{14mu} 8} \right\rbrack\end{matrix}$ where N: a total number of pixels in the local regionimage, N_(i): a total number of pixels whose pixel value is a digitalvalue of i, N_(j): a total number of pixels whose pixel value is adigital value of j k: a constant, a pixel value is changed from thedigital value i to the digital value j by the modification by the means(b).
 24. The image reconstructing device of claim 20, wherein the means(h) calculates ΔS using [Expression 7], and obtains ΔE (ΔE=ΔH−TΔS+ . . .) including as a component a sum of a product (−TΔS) of the calculatedΔS and the temperature parameter (T), and the ΔH,ΔS=k ln N _(i) −k ln(N _(j)+1)  [Expression 7] where S represents anentropy of a local region image of d×d pixels around the coordinatepoint (x₀, y₀) and is calculated by [Expression 8], $\begin{matrix}{S = {k\; \ln \frac{N!}{{N_{1}!}{N_{2}!}\mspace{14mu} \ldots \mspace{14mu} {N_{i}!}\mspace{20mu} \ldots \mspace{14mu} {N_{n}!}}}} & \left\lbrack {{Expression}\mspace{14mu} 8} \right\rbrack\end{matrix}$ where N: a total number of pixels in the local regionimage, N_(i): a total number of pixels whose pixel value is a digitalvalue of i, N_(j): a total number of pixels whose pixel value is adigital value of j, k: a constant, a pixel value is changed from thedigital value i to the digital value j by the modification by the means(b).
 25. The image reconstructing device of claim 15, comprising:instead of the means (e) and (f), means (e1) for determining whether ornot the modification is to be accepted, based on an acceptance functionusing the differential (ΔE) and a temperature parameter (T) forcontrolling an acceptance probability, and reserving reflection of aresult of determination on the current estimated cross-sectional image;and means (f1) for reflecting the reservation(s) in the means (e1) onthe current estimated cross-sectional image and changing a value of thetemperature parameter (T) every time the number of iterations of aseries of processes of the means (a) to (d) and (e1) reaches apredetermined value.
 26. An image reconstructing device for obtaining across-sectional image of an object from projections obtained byirradiating the object with a beam of radiation, comprising: means (m1)for calculating a back projection g₀(x, y) of a radiographic projectionp₀(r, θ) of the object by a back projection operation without filtering;means (m2) for calculating a projection p(r, θ) from a current estimatedcross-sectional image f(x, y) of the object, and calculating a backprojection g(x, y) of the projection p(r, θ) by a back projectionoperation without filtering; means (m3) for generating an image Δμ(x, y)whose pixel value is a change value of the current estimatedcross-sectional image f(x, y) of the object; means (m4) for generatingan image ΔH(x, y) by applying [Expression 9] to each pixel value,ΔH=MΔμ ²+2Δμ{g(x ₀ ,y ₀)−g ₀(x ₀ ,y ₀)}  [Expression 9] where M: thenumber of projection angles; means (m5) for calculating ΔE(x, y) usingthe ΔH(x, y), where ΔE(x, y) represents a differential betweenevaluation functions E₀(x, y) and E₁(x, y), E₀(x, y) represents anevaluation function including a difference between the projection p(r,θ) calculated from the estimated cross-sectional image f(x, y) and theradiographic projection p₀(r, θ), and E₁(x, y) represents an evaluationfunction including a difference between a projection {p(r, θ)+Δp(r, θ)},calculated from a sum {f(x, y)+Δμ(x, y)} of the estimatedcross-sectional image f(x, y) and the image Δμ(x, y) obtained by themeans (m3), and the radiographic projection p₀(r, θ); means (m6) forsetting the Δμ(x, y) to 0 at a coordinate point (x, y) where the ΔE ispositive; and means (m7) for setting a sum of the estimatecross-sectional image f(x, y) and the image Δμ(x, y) obtained by themeans (m6) as a new estimated cross-sectional image f(x, y) andrepeating processes of the means (m2) to (m6) with respect to the newestimated cross-sectional image f(x, y).
 27. An image reconstructingdevice for obtaining a cross-sectional image of an object fromprojections obtained by irradiating the object with a beam of radiation,comprising: means (m1) for calculating a back projection g₀(x, y) of aradiographic projection p₀(r, θ) of the object using [Expression 10]$\begin{matrix}{{g_{0}\left( {x,y} \right)} = {\sum\limits_{\theta}{p_{0}\left( {{{x\; \cos \; \theta} + {y\; \sin \; \theta}},\theta} \right)}}} & \left\lbrack {{Expression}\mspace{14mu} 10} \right\rbrack\end{matrix}$ r: a channel position of a one-dimensional detector takingthe projection, θ: a projection angle, means (m2) for calculating aprojection p(r, θ) from a current estimated cross-sectional image f(x,y) of the object, and calculating a back projection g(x, y) of theprojection p(r, θ) using [Expression 11]${g\left( {x,y} \right)} = {\sum\limits_{\theta}{p\left( {{{x\; \cos \; \theta} + {y\; \sin \; \theta}},\theta} \right)}}$means (m3) for generating an image Δμ(x, y) whose pixel value is achange value of the current estimated cross-sectional image f(x, y) ofthe object; means (m4) for generating an image ΔH(x, y) by applying[Expression 12] to each pixel value,ΔH=MΔμ ²+2Δμ{g(x ₀ ,y ₀)−g ₀(x ₀ ,y ₀)}  [Expression 12] where M: thenumber of projection angles; means (m5) for generating Δσ(x, y) byapplying [Expression 13] to each pixel value, $\begin{matrix}{{\Delta \; \sigma} = \sqrt{\begin{matrix}{{\langle{f\left( {x_{0},y_{0}} \right)^{2}}\rangle} - \frac{f_{i}^{2} - f_{j}^{2}}{d^{2}} -} \\\left\{ {{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} - \frac{f_{i} - f_{j}}{d^{2}}} \right\}^{2}\end{matrix}}} & \left\lbrack {{Expression}\mspace{14mu} 13} \right\rbrack\end{matrix}$ where σ represents a standard deviation of luminancevalues of d×d pixels around a coordinate point (x₀, y₀) and iscalculated by [Expression 14], and f_(i) and f_(j) represent values off(x₀, y₀) before and after a change, $\begin{matrix}{{\sigma = \sqrt{{\langle{f\left( {x_{0},y_{0}} \right)}^{2}\rangle} - {\langle{f\left( {x_{0},y_{0}} \right)}\rangle}^{2}}}{where}} & \left\lbrack {{Expression}\mspace{14mu} 14} \right\rbrack \\{{\langle{f\left( {x_{0},y_{0}} \right)}^{2}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{f\left( {{x_{0} + i},{y_{0} + j}} \right)}^{2}}}} & \left\lbrack {{Expression}\mspace{14mu} 15} \right\rbrack \\{{\langle{f\left( {x_{0},y_{0}} \right)}\rangle} = {\frac{1}{d^{2}}{\sum\limits_{i,{j = {{- d}/2}}}^{d/2}{f\left( {{x_{0} + i},{y_{0} + j}} \right)}}}} & \left\lbrack {{Expression}\mspace{14mu} 16} \right\rbrack\end{matrix}$ means (m6) for generating an image ΔS by applying[Expression 17] to each value,ΔS=k ln N _(i) −k ln(N _(j)+1)  [Expression 17] where S represents anentropy of a local region image of d×d pixels around the coordinatepoint (x₀, y₀) and is calculated by [Expression 18], $\begin{matrix}{S = {k\; \ln \frac{N!}{{N_{1}!}{N_{2}!}\mspace{14mu} \ldots \mspace{14mu} {N_{i}!}\mspace{20mu} \ldots \mspace{14mu} {N_{n}!}}}} & \left\lbrack {{Expression}\mspace{14mu} 18} \right\rbrack\end{matrix}$ where N: a total number of pixels in the local regionimage, N_(i): a total number of pixels whose pixel value is a digitalvalue of i, N_(j): a total number of pixels whose pixel value is adigital value of j, k: a constant, a pixel value is changed from thedigital value i to the digital value j by the modification by the means(b); means (m7) for calculating ΔE(x, y) based on [Expression 19],ΔE=ΔH+cΔσ−TΔS  [Expression 19] where c: a coefficient, T: a virtualtemperature (temperature parameter), means (m8) for setting the Δt(x, y)to 0 at a coordinate point (x, y) where the ΔE is positive; means (m9)for setting a sum of the estimated cross-sectional image f(x, y) and theimage Δμ(x, y) obtained by the means (m8) as a new estimatedcross-sectional image f(x, y); and means (m10) for multiplying the T byα (α<1), and repeating processes of the means (m2) to (m9).
 28. A methodfor obtaining a cross-sectional image of an object from projections(hereinafter referred to as “radiographic projections”) obtained byirradiating the object with a beam of radiation, comprising the stepsof: (a) obtaining an evaluation function (hereinafter referred to as an“energy”) (E₀) including differences between projections calculated froma current estimated cross-sectional image of the object and theradiographic projections; (b) modifying a portion of the currentestimated cross-sectional image; (c) obtaining an energy (E₁) includingdifferences between projections calculated from the modified estimatedcross-sectional image and the radiographic projections; (d) obtaining adifferential (ΔE) between the energy (E₀) and the energy (E₁); (e)determining whether or not the modification is to be accepted, based onan acceptance function using the differential (ΔE) and a temperatureparameter (T) for controlling an acceptance probability; (f) reflectinga result of the determination on the current estimated cross-sectionalimage, and returning to the step (a); (g) changing a value of thetemperature parameter (T) every time the number of iterations of thesteps (a) to (f) reaches a predetermined value; and (h) determiningwhether or not the result of the determination in the step (e) satisfiespredetermined stop conditions, and if the result of the determination inthe step (e) satisfies predetermined stop conditions, ending theprocess.
 29. An image reconstructing program for obtaining across-sectional image of an object from projections (hereinafterreferred to as “radiographic projections”) obtained by irradiating theobject with a beam of radiation, wherein the program causes a computerto execute the steps of: (a) obtaining an evaluation function(hereinafter referred to as an “energy”) (E₀) including differencesbetween projections calculated from a current estimated cross-sectionalimage of the object and the radiographic projections; (b) modifying aportion of the current estimated cross-sectional image; (c) obtaining anenergy (E₁) including differences between projections calculated fromthe modified estimated cross-sectional image and the radiographicprojections; (d) obtaining a differential (ΔE) between the energy (E₀)and the energy (E₁); (e) determining whether or not the modification isto be accepted, based on an acceptance function using the differential(ΔE) and a temperature parameter (T) for controlling an acceptanceprobability; (f) reflecting a result of the determination on the currentestimated cross-sectional image, and returning to the step (a); and (g)changing a value of the temperature parameter (T) every time the numberof iterations of the steps (a) to (f) reaches a predetermined value. 30.A CT apparatus comprising: means (A) for obtaining projections byirradiating an object with a beam of radiation; and means (B) forobtaining a cross-sectional image of the object from the projections,wherein the means (B) includes: means (b1) for obtaining an evaluationfunction (hereinafter referred to as an “energy”) (E₀) includingdifferences between projections calculated from a current estimatedcross-sectional image of the object and the projections by irradiatingthe object with the beam of radiation (hereinafter referred to as“radiographic projections”); means (b2) for modifying a portion of thecurrent estimated cross-sectional image; means (b3) for obtaining anenergy (E₁) including differences between projections calculated fromthe modified estimated cross-sectional image and the radiographicprojections; means (b4) for obtaining a differential (ΔE) between theenergy (E₀) and the energy (E₁); means (b5) for determining whether ornot the modification is to be accepted, based on an acceptance functionusing the differential (ΔE) and a temperature parameter (T) forcontrolling an acceptance probability, and reflecting a result of thedetermination on the current estimated cross-sectional image; and means(b6) for changing a value of the temperature parameter (T) every timethe number of iterations of a series of processes of the means (b1) to(b5) reaches a predetermined value.